// 版权所有2009 Go作者。版权所有。
// 此源代码的使用受BSD样式
// 许可证的约束，该许可证可以在许可证文件中找到。

package math

func isOddInt(x float64) bool {
	xi, xf := Modf(x)
	return xf == 0 && int64(xi)&1 == 1
}

// 特殊情况取自FreeBSD的/usr/src/lib/msun/src/e_pow.c 
// 由IEEE标准754-2008“第9.2.1节特殊值”更新。

// Pow返回x**y，y的基x指数。
// 
// 特殊情况是（按顺序）：
// Pow（x，±0）=1表示任意x 
// Pow（1，y）=1表示任意y 
// Pow（x，1）=x表示任意x 
// Pow（NaN y）=NaN NaN 
// Pow（x，NaN NaN NaN NaN 
// Pow（±0，y）=±Inf表示y为奇整数<0 
// Pow（±0，-Inf）=+Inf 
// Pow（±0，+Inf）=+0 
// Pow（±0，y）=+Inf表示有限y<0而非奇整数
// Pow（±0，y）=±0表示y为奇整数>0 
// Pow（±0，y）=+0表示有限y>0，而不是奇数整数
// /Pow（-1，±Inf）=1 
// /Pow（x，+Inf）=+Inf表示| x |>1 
// Pow（+Inf，y）=+0表示y<0 
// Pow（-Inf，y）=Pow（-0，-y）
// Pow（x，y）=NaN表示有限x<0和有限非整数y 
func Pow(x, y float64) float64 {
	if haveArchPow {
		return archPow(x, y)
	}
	return pow(x, y)
}

func pow(x, y float64) float64 {
	switch {
	case y == 0 || x == 1:
		return 1
	case y == 1:
		return x
	case IsNaN(x) || IsNaN(y):
		return NaN()
	case x == 0:
		switch {
		case y < 0:
			if isOddInt(y) {
				return Copysign(Inf(1), x)
			}
			return Inf(1)
		case y > 0:
			if isOddInt(y) {
				return x
			}
			return 0
		}
	case IsInf(y, 0):
		switch {
		case x == -1:
			return 1
		case (Abs(x) < 1) == IsInf(y, 1):
			return 0
		default:
			return Inf(1)
		}
	case IsInf(x, 0):
		if IsInf(x, -1) {
			return Pow(1/x, -y) // Pow（-0，-y）
		}
		switch {
		case y < 0:
			return 0
		case y > 0:
			return Inf(1)
		}
	case y == 0.5:
		return Sqrt(x)
	case y == -0.5:
		return 1 / Sqrt(x)
	}

	yi, yf := Modf(Abs(y))
	if yf != 0 && x < 0 {
		return NaN()
	}
	if yi >= 1<<63 {
		// yi是一个大偶数，将导致溢出（或下溢到0）
		// 对于除-1之外的所有x（x==1在前面处理过）
		switch {
		case x == -1:
			return 1
		case (Abs(x) < 1) == (y > 0):
			return 0
		default:
			return Inf(1)
		}
	}

	// ans=a1*2**ae（=现在为1）。
	a1 := 1.0
	ae := 0

	// ans*=x**yf 
	if yf != 0 {
		if yf > 0.5 {
			yf--
			yi++
		}
		a1 = Exp(yf * Log(x))
	}

	// ans*=x**yi 
	// 根据yi位乘以x的连续平方
	// 。
	// 将二的幂累加到exp.
	x1, xe := Frexp(x)
	for i := int64(yi); i != 0; i >>= 1 {
		if xe < -1<<12 || 1<<12 < xe {
			// 在xe溢出左移位之前捕获xe，
			// 因为我=0它至少还设置了一个位，因此ae将在至少一次迭代中累积xe 
			// ae+=xe是ae 
			// 的下限，ae的下限超过浮点64 exp 
			// 因此对Ldexp的最终调用将产生under/overflow（0/Inf）
			ae += xe
			break
		}
		if i&1 == 1 {
			a1 *= x1
			ae += xe
		}
		x1 *= x1
		xe <<= 1
		if x1 < .5 {
			x1 += x1
			xe--
		}
	}

	// ans=a1*2**ae 
	// 如果y<0{ans=1/ans}
	// 但顺序相反
	if y < 0 {
		a1 = 1 / a1
		ae = -ae
	}
	return Ldexp(a1, ae)
}
